* Define program to estimate the damages
* This program is used inside the next program "ricardian_program"
capture program drop damages
program define damages, rclass

	syntax [, deficit(string)  surplus(string) edd(string) ///
			knots_d(real 0) knots_s(real 0) knots_g(real 0) knots_e(real 0) ///
			practice(real 0) rent_hat(string)]
	qui {
	local knots_`deficit' `knots_d'
	local knots_`surplus' `knots_s'
	local knots_gdd `knots_g'
	local knots_`edd' `knots_e'

	* Use same knot locations as in my regression
	foreach var of varlist `deficit' `surplus' gdd `edd' {
	local knot1=`var'_knots[1,1]
	local knot2=`var'_knots[1,2]
	local knot3=`var'_knots[1,3]
	local knot4=`var'_knots[1,4]
	local knot5=`var'_knots[1,5]
	local knot6=`var'_knots[1,6]
	capture drop `var'_spline*
		if `knots_`var''==2 {
		gen `var'_spline = `var'
		}
		if `knots_`var''==3 {
		mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3')
		}
		if `knots_`var''==4 {
		mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3'  `knot4')
		}
		if `knots_`var''==5 {
		mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3'  `knot4' `knot5')
		}
		if `knots_`var''==6 {
		mkspline `var'_spline = `var', cubic knots(`knot1'  `knot2'  `knot3'  `knot4' `knot5' `knot6')
		}
	}
	
	estimates restore main_results
	if `practice'==0 {
	predict ln_rent_hat_new, xb
	}
	if `practice'==1 {
	predict ln_rent_hat_new if ln_cash_rent!=., xb
	}
	gen rent_hat_new=exp(ln_rent_hat_new+e(rmse)^2/2)
	
	gen chg_rent=rent_hat_new - `rent_hat'
	gen ln_chg_rent= ln_rent_hat_new - ln_rent_hat
	gen perc_chg_rent=chg_rent/`rent_hat'
	
	* Total Damages
	if `practice'==0 {
		gen damages=chg_rent*croplandNon_acres
	}
	if `practice'==1 {
		* make sure change estimates are missing if no irrigated acres harvested
		replace chg_rent=. if acresh_Irr2012==0 | acresh_Irr2012==.
		replace ln_chg_rent=. if acresh_Irr2012==0 | acresh_Irr2012==.
		replace perc_chg_rent=. if acresh_Irr2012==0 | acresh_Irr2012==.
		gen damages=chg_rent*acresh_Irr2012
	}
	total damages
	return scalar total_damages=_b[damages]
	} // end of qui

end


capture program drop ricardian_program
program ricardian_program, eclass
tempname parameters 

	syntax [, deficit(string)  surplus(string) edd(string) ///
			knots_d(real 0) knots_s(real 0) knots_g(real 0) knots_e(real 0) ///
			practice(real 0) reg_year(real 0) soils_controls(string) scenario(string) model(string) wild(real 0)]

tempvar rent_hat total_rent
local main_year 2013

if `wild'==1{
capture drop temp pos wildresid wildy
bys stfips_clust: gen temp = uniform() 
by stfips_clust: gen pos = (temp[1] < .5) 
gen wildresid = epshat * (2*pos - 1) 
gen wildy = xb + wildresid 
}		

if `knots_d'==2 {
capture drop `deficit'_spline*
gen `deficit'_spline = `deficit'
matrix `deficit'_knots=[0]
}

if `knots_d'>=3 {
capture drop `deficit'_spline*
mkspline `deficit'_spline = `deficit', cubic nknots(`knots_d') 
matrix `deficit'_knots=r(knots)
}

if `knots_s'==2 {
capture drop water_surplus_spline*
gen `surplus'_spline = `surplus'
matrix `surplus'_knots=[0]
}

if `knots_s'>=3 {
capture drop water_surplus_spline*
mkspline `surplus'_spline = `surplus', cubic nknots(`knots_s') 
matrix `surplus'_knots=r(knots)
}

if `knots_g'==2 {
capture drop gdd_spline*
gen gdd_spline = gdd
matrix gdd_knots=[0]
}

if `knots_g'>=3 {
capture drop gdd_spline*
mkspline gdd_spline = gdd, cubic nknots(`knots_g') 
matrix gdd_knots=r(knots)
}

if `knots_e'==2 {
capture drop `edd'_spline*
gen `edd'_spline = `edd'
matrix `edd'_knots=[0]	
}

if `knots_e'>=3 {
capture drop `edd'_spline*
mkspline `edd'_spline = `edd', cubic nknots(`knots_e') 
matrix `edd'_knots=r(knots)		
}

if `wild'==0 & `practice'==0 {			
reg ln_cash_rent `deficit'_spline* `surplus'_spline* gdd_spline* `edd'_spline* `soils_controls' i.stfips ///
			[aw=croplandNon_acres], vce(cluster stfips_clust)
est sto main_results
* Predicted value and residual used for wild bootstrap
qui predict xb if e(sample), xb
qui predict epshat if e(sample), resid
}

if `wild'==0 & `practice'==1 {			
reg ln_cash_rent `deficit'_spline* `surplus'_spline* gdd_spline* `edd'_spline* `soils_controls' i.stfips ///
			[aw=acresh_Irr2012], vce(cluster stfips_clust)
est sto main_results
* Predicted value and residual used for wild bootstrap
qui predict xb if e(sample), xb
qui predict epshat if e(sample), resid
}

if `wild'==1 & `practice'==0 {			
reg wildy `deficit'_spline* `surplus'_spline* gdd_spline* `edd'_spline* `soils_controls' i.stfips ///
			[aw=croplandNon_acres], vce(cluster stfips_clust)
est sto main_results
}

if `wild'==1 & `practice'==1 {			
reg wildy `deficit'_spline* `surplus'_spline* gdd_spline* `edd'_spline* `soils_controls' i.stfips ///
			[aw=acresh_Irr2012], vce(cluster stfips_clust)
est sto main_results
}

qui {
*-------------------------------------------------------------------------------
* Estimate Climate Change Damages
*-------------------------------------------------------------------------------
if `practice'==0{
predict ln_rent_hat, xb
}
* Only predict irrigated rent if had irrigated cash rent data
if `practice'==1{
predict ln_rent_hat if ln_cash_rent!=., xb
}
gen `rent_hat'=exp(ln_rent_hat+e(rmse)^2/2)
drop ln_rent_hat
* total rent in the region
if `practice'==0 {
	gen `total_rent'=`rent_hat'*croplandNon_acres
}
if `practice'==1 {
	gen `total_rent'=`rent_hat'*acresh_Irr2012
}

total `total_rent'
scalar base_rent=_b[`total_rent']
scalar ln_base_rent=ln(base_rent)


*******************************
* Total Damages
*******************************
preserve
	foreach var of varlist `deficit' `surplus' gdd `edd' {
		summ `var' chg_`var'_rcp45AVG
		replace `var'=`var' + chg_`var'_`scenario'`model'
	}
	replace water_deficitOffSeason=water_deficitOffSeason + chg_deficitOff_`scenario'`model'
	replace water_surplusOffSeason=water_surplusOffSeason + chg_surplusOff_`scenario'`model'
	damages, deficit("`deficit'") surplus("`surplus'") edd("`edd'") ///
			knots_d(`knots_d') knots_s(`knots_s') knots_g(`knots_g') knots_e(`knots_e') ///
			practice(`practice') rent_hat(`rent_hat')
	return list, all
	* Save damage estimates for map of damages
	* Only save if the main estimation and not bootstrap estimation
	if `wild'==0 & `reg_year'==`main_year' & `practice'==0 {	
	save "..\temp\damages_total_`scenario'`model'", replace
	}
	if `wild'==0 & `reg_year'==`main_year' & `practice'==1 {	
	save "..\temp\damages_irr_total_`scenario'`model'", replace
	}
restore

* Relative change in rent
local perc_chg_rent_total=r(total_damages)/base_rent
display "`perc_chg_rent_total'"
local damages_total=r(total_damages)
* Change in log rent
local ln_new_rent=r(total_damages)+base_rent
local chg_ln_rent_total= ln(`ln_new_rent') - ln_base_rent
display "`chg_ln_rent_total'"
* Relative change in rent implied by log change
display exp(`chg_ln_rent_total') - 1

*------------------------------------------------------------------------------------
/* If I calculate the percent change in rent due to climate damages
   for each variable separately, then the perecent changes do not all
   add to 1 perfectly due to the log tranformation. But note that the change
   in natural log only approximates a percent change and this approximation
   is less accurate for larger changes. Above I have shown two alternative 
   ways of calculating the percent change in rent and get the same answer either way.
   To decompose the total change in rent and find the share of damages due to each
   source, I calculate the share of the total change in log rent because these 
   shares add to 1. */
*------------------------------------------------------------------------------------

*******************************
* Only Water Deficit 
*******************************
preserve
	foreach var of varlist `deficit'  {
		replace `var'=`var' + chg_`var'_`scenario'`model'
	}
	replace water_deficitOffSeason=water_deficitOffSeason + chg_deficitOff_`scenario'`model'
	damages, deficit("`deficit'") surplus("`surplus'") edd("`edd'") ///
			knots_d(`knots_d') knots_s(`knots_s') knots_g(`knots_g') knots_e(`knots_e') ///
			practice(`practice') rent_hat(`rent_hat')
	return list, all
	if `wild'==0 & `reg_year'==`main_year' & `practice'==0 {	
	save "..\temp\damages_deficit_`scenario'`model'", replace
	}
	if `wild'==0 & `reg_year'==`main_year' & `practice'==1 {	
	save "..\temp\damages_irr_deficit_`scenario'`model'", replace
	}
restore

* Relative change in rent
local perc_chg_rent_deficit=r(total_damages)/base_rent
local damages_deficit=r(total_damages)
* Proportion due to deficit
local share_deficit = (ln(r(total_damages)+base_rent) - ln_base_rent) /`chg_ln_rent_total'

*******************************
* Only Water Surplus
*******************************
preserve
	foreach var of varlist `surplus'  {
		replace `var'=`var' + chg_`var'_`scenario'`model'
	}
	replace water_surplusOffSeason=water_surplusOffSeason + chg_surplusOff_`scenario'`model'
	damages, deficit("`deficit'") surplus("`surplus'") edd("`edd'") ///
			knots_d(`knots_d') knots_s(`knots_s') knots_g(`knots_g') knots_e(`knots_e') ///
			practice(`practice') rent_hat(`rent_hat')
	return list, all
	if `wild'==0 & `reg_year'==`main_year' & `practice'==0 {	
	save "..\temp\damages_surplus_`scenario'`model'", replace
	}
	if `wild'==0 & `reg_year'==`main_year' & `practice'==1 {	
	save "..\temp\damages_irr_surplus_`scenario'`model'", replace
	}
restore

* Relative change in rent
local perc_chg_rent_surplus=r(total_damages)/base_rent
local damages_surplus=r(total_damages)
* Proportion due to surplus
local share_surplus = (ln(r(total_damages)+base_rent) - ln_base_rent) /`chg_ln_rent_total'

*******************************
* Only Heat Stress
*******************************
preserve
	foreach var of varlist  gdd `edd' {
		replace `var'=`var' + chg_`var'_`scenario'`model'
	}
	damages, deficit("`deficit'") surplus("`surplus'") edd("`edd'") ///
			knots_d(`knots_d') knots_s(`knots_s') knots_g(`knots_g') knots_e(`knots_e') ///
			practice(`practice') rent_hat(`rent_hat')
	return list, all
	if `wild'==0 & `reg_year'==`main_year' & `practice'==0 {	
	save "..\temp\damages_heat_`scenario'`model'", replace
	}
	if `wild'==0 & `reg_year'==`main_year' & `practice'==1 {	
	save "..\temp\damages_irr_heat_`scenario'`model'", replace
	}
restore

* Relative change in rent
local perc_chg_rent_heat=r(total_damages)/base_rent
local damages_heat=r(total_damages)
* Proportion due to heat
local share_heat = (ln(r(total_damages)+base_rent) - ln_base_rent) /`chg_ln_rent_total'



} // end of qui

display "Total predicted damages " `damages_total'
display "Water Deficit predicted damages " `damages_deficit'
display "Water Surplus predicted damages " `damages_surplus'
display "Heat Stress predicted damages " `damages_heat'

matrix `parameters'= (`perc_chg_rent_total',`perc_chg_rent_deficit',`perc_chg_rent_surplus',`perc_chg_rent_heat', ///
					`share_deficit',`share_surplus',`share_heat')

matrix colnames `parameters' = rel_total_chg rel_deficit_chg rel_surplus_chg rel_heat_chg ///
				share_deficit share_surplus share_heat 

ereturn post `parameters'
end
